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Abstract 

We investigate Isotropic - Nematic transition in liquid crystal elastomers employing non- 
Boltzmann Monte Carlo techniques. We consider a lattice model of a liquid elastomer and a recently 
proposed Hamiltonian which accounts for homogeneous/inhomogeneous interactions among liquid 
crystalline units, interaction of local nematics with global strain, and with inhomogeneous fields 
and stress. We find that when the local director is coupled strongly to the global strain the tran- 
sition is strongly first order; the transition becomes weakly first order when the coupling becomes 
weaker. The transition temperature decreases with decrease of coupling strength. Besides, we find 
that the nematic order scales nonlinearly with global strain especially for strong coupling and at 
low temperatures. 

PACS numbers: 64.70.Md.61.30.Vx; 61.30.Cz; 61.41. +c 
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I. INTRODUCTION 



A liquid crystal elastomer [l| is a weakly cross-linked, percolating network of long, rigid, 
liquid crystalline units. It inherits and combines the properties of both of its components - 
rubber elasticity of the cross-linked network and nematic and smectic ordering capabilities of 
the liquid crystalline units. As a result, a liquid crystal elastomer exhibits unusual properties 

n h, nnnnnu 

that are of interest from both basic |2J, y( and applied science |4|, [5|, |6|, |7|, |8|, 12[ points of views. 
For example even a small change in temperature (that causes the liquid crystalline units 
to transform from nematic to isotropic phase) can lead to large and reversible changes in 
the shape of a liquid crystal elastomer 3|, |4( . It is precisely this property that makes it a 
competitive candidate material for artificial muscles, see e.g. 

aa. 

In fact, since liquid 

crystal elastomers respond sensitively, not only to temperature, but also to other stimuli, 
like electric fields, magnetic fields, ultra violet radiations, gamma rays etc., they become 
suitable for construction of actuators, detectors, micro-pumps etc., see e.g. 0,0,3,3. 

We have carried out Monte Carlo simulation of a lattice model, employing non-Boltzmann 
sampling techniques. The transition from disordered isotropic phase at high temperatures 
to ordered nematic phase at low temperatures is known to be weakly first order, in bulk 
liquid crystalline material, see e.g. [10|]. Our simulations show that in a liquid crystal 
elastomer, when the local nematic is coupled strongly to global elastic order the transition 
is strongly first order; the transition becomes weak when the coupling strength decreases. 
It is known that the above model exhibits strong discontinuity for 7 = 1 [11] and a weak 
discontinuity for 7 = 0, see 12] . Our work interpolates these two limiting cases. More 
importantly, we find that the nematic order S, scales non-linearly with the strain parameter 
A, for strong coupling and at low temperatures. Our Monte Carlo simulations are based on a 



modified 



13] Wang-Landau implementation [l^] of entropic sampling techniques [15|] useful 



for studying complex systems in the neighbourhood of phase transition. In general, entropic 
sampling 15[ and related techniques 14J enable calculation of the macroscopic properties 
of a closed system over a wide range and at arbitrarily fine resolution of temperatures - all 
from a single entropic ensemble that contains microstates of all energies in approximately 
equal proportions. More importantly these techniques help estimate free energy and entropy, 
a task which is very difficult, if not impossible, in conventional Metropolis Markov chain 
Monte Carlo techniques. From the plots of free energy versus order parameter at various 
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temperatures, we can unambiguously determine the nature of the phase transition. The 
paper is organized as follows. In section II we describe a lattice model of liquid crystal 



elastomer and the Hamiltonian, proposed by Selinger and Ratna 11]. In section III we 
describe the non-Boltzmann Monte Carlo simulation techniques employed for studying the 
problem. Section IV is devoted to results and discussions. Principal outcomes of the study 
are brought out briefly in the concluding section V. 

II. LATTICE MODEL OF A LIQUID CRYSTAL ELASTOMER 

attice model of a liquid crystal elastomer system proposed by Selinger, 



We employ a 

Jeon and Ratna Uj. We consider an L x L x L cubic lattice with lattice sites indexed 
by natural numbers i. Each lattice site holds a nematic director denoted by a unit vector 
\ui). Each nematic director interacts with its six nearest neighbours and the interaction 
is described by Lebwohl-Lasher potential [161 ]. which has head-tail flip symmetry (\ui) and 
—\ui) are equivalent). Besides, each director is coupled to global elastic degree of freedom. 
A model Hamiltonian for such an interaction between local director and global strain has 
been derived starting from the neo-classical theory of rubber elasticity [l|, see ll|. The 



Hamiltonian includes 

(i) Lebwohl-Lasher nearest neighbour interaction [16] of the directors, placed on the lattice 
sites, 

(ii) the interaction of each local director with 

(a) the global elastic degree of freedom and 

(b) an inhomogeneous field, 

(iii) an externally imposed global stress and 

(iv) shear modulus that couples to the strain and to local directors. 
It is given by, 
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E = - Yl Ji >i ( \ ( u i\ u i) 2 



L 3 

£ 

i=i 



X 2 V A / 2 2 



In the above the first term is the Lebwohl - Lasher potential. The sum extends over all 
distinct pairs (i, j) of nearest neighbours. Jij > measures the strength of nearest neighbour 
interaction; periodic boundary condition is imposed in all directions, /i is the shear modulus. 
For a given configuration (microstate C) of the directors, we first construct an average 
projection operator given by 

L 3 

A(C) = ^^ i \u i )(u i \ • ( 2 ) 
i=i 

From A(C) we construct a traceless symmetric tensor, 

Q(C) = A(G) - ^trace(A) x I , (3) 

where / denotes a unit matrix. Let rj denote the largest eigenvalue of Q. The orientational 
order parameter of the system when in microstate C is given by S = Srj/2. The corresponding 
eigenvector is denoted by \m). X > 1 is a scalar denoting the strain parameter. We take 
A = 1 + e, where e is the strain along the distortion axis taken to be along the vector \m). 
The applied stress is denoted by a and \hi) denotes an inhomogeneous field experienced by 
the director at lattice site i. The coupling of the nematic to the elastic degree of freedom is 
tuned by the parameter < 7 < 1. The above model Hamiltonian can describe liquid crystal 
elastomer with random bond {{Jij}) and random field disorder, in the presence of 

external stress, o 7^ 0. 



Selinger and Ratna [ll| considered first a homogeneous elastomer (Jjj = 1 V ) in the 
absence of local fields (\hi) — V i) with maximal coupling (7 = 1) and no external stress 
(a = 0). We observe that setting the coupling 7 = 1, especially at low temperatures, is not 
realistic. For, such a choice would correspond to a very steep anisotropy and would imply 
extreme elongation in one direction. On the other hand if 7 = 0, the nematic order would 
get decoupled from the global strain. Phase transitions in the liquid crystalline units will 
have no effect on the shape of the elastomer. The actual scenario is likely to correspond to 



4 



a choice of 7 between zero and unity. Accordingly, in our study we consider a homogeneous 
system (Jjj = 1 V with no external fields = V i) and without any applied 

stress (a = 0). We set // = 1 as recommended in [111]. This corresponds to setting the 
network energy at zero strain to the highest possible energy of the liquid crystalline units. 
The network is robust and deforms only under large stress. 

III. MONTE CARLO SIMULATION 

A typical Markov chain Monte Carlo simulation of the system would proceed as follows. 
Start with an initial microstate Co of the system defined by a configuration of directors 
and strain parameter A > 1. Let Eq denote the energy of the microstate, calculated 
from Eq. ([!]). Select randomly a director; rotate it randomly so that it points in a direction 



within a specified cone about its current direction. Barker's method [17| for selecting the 
trial orientation is usually employed. Change the value of A randomly. We thus have a trial 
microstate Ct of the liquid crystal elastomer system. Let E t be its energy. Accept the trial 
state with a probability 



p = min ^1, exp [ — (3 AE ] J . 



(4) 



where AE = E t — Eq. Here (3 = 1/[/cbT] with T denoting temperature and ks, the Boltz- 
mann constant, set to unity. Thus, the next microstate Ci is either C t with probability p 
or Co with probability 1 — p. Proceed in the same fashion and construct a Markov chain of 
microstates denoted by Co — > Ci — > C2 — > • • • C n — > ■ ■ ■ . This is called the Metropolis algo- 
rithm The asymptotic (n — > 00) part of the Markov chain would contain microstates 
that belong to a canonical ensemble at the temperature chosen for the simulation. 

Metropolis algorithm and it variants, see e.g. 19(, come under the class of Boltzmann 
sampling techniques. The limitations of Boltzmann sampling have long since been recog- 
nized. For example it can not address satisfactorily problems of super-critical slowing down 
near first order phase transitions, an issue of relevance to the simulations of liquid crystal 
elastomer systems considered in this paper. The microstates representing the interface be- 
tween the ordered and disordered phases have intrinsically low probability of occurrence in 
a closed system and hence are scarcely sampled. Switching from one phase to the other 
takes a very long time due to presence of high energy barriers when the system size is large. 



5 



As a result the relative free energies of ordered and disordered phases can not be easily and 
accurately determined. 

For these kinds of problems, non-Boltzmann sampling provides a legitimate alternative. 
Accordingly, for the simulation of lattice elastomers, we consider entropic sampling 15], 
which is based on the following premise. The probability that a closed system can be found 
in a microstate C is given by 

P(C) = [Z{P)]- l eM-PE{C)} (5) 

where E{C) is the energy of the microstate C, and Z{(5) is the canonical partition function 
given by, 

Z{fi) = J>xp[-/3£(C)] 
c 

= J2 D ( E ) e M-PE}. (6) 

E 

In the above D(E) is the density of states. The probability density for a closed system to 
have an energy E is thus proportional to D(E) exp[—@E]. Let us suppose we want to sample 
microstates in such a way that the resultant probability density of energy is 

P 9 (E) cc D(E) [g(E)]- 1 (7) 

where g(E) > V E, is a function of your choice. Non-Boltzmann sampling of microstates, 
consistent with Eq. (j7j) is implemented as follows. Let Q be the current and C t the trial 
microstates respectively. Let E{ = E(Ci) and E t = E(C t ) denote the energy of the current 
and of the trial microstates respectively. The next entry Gj + i in the Markov chain is taken 
as, Ct with probability p and Q with probability 1 — p, and p is given by, 

g{Ei 



p = mm 



(8) 



g(E t 

It is easily verified that the above acceptance rule obeys detailed balance and hence we are 
assured that the Markov chain constructed would converge asymptotically to the desired 
g ensemble. When [g(E)]^ 1 = exp(—/3E) we recover conventional Boltzmann sampling 
implemented in the Metropolis algorithm, described earlier. For any other choice of g(E) 
we get non-Boltzmann sampling. However, canonical ensemble average of a macroscopic 
property 0(C) can be obtained by un- weighting and re- weighting of 0(C) for each C sampled 



from the g ensemble; for un-weighting we divide by [g(E(C))] 1 and for re-weighting we 
multiply by exp[—/3E(C)}. Formally we have, 

(n) E C 0(C)g(E(C))exp[-(3E(C)] . , 

{ } Ec9(E(C))ex P [-(3E(C)] ' 1 ) 

The left hand side of the above is the equilibrium value of O in a closed system at (3, 
while in the right side, the summation in the numerator and in the denominator, runs over 
microstates belonging to the non-Boltzmann g ensemble. The important point is that if 
we take a suitable and temperature independent g, then from a single g ensemble we can 
calculate the macroscopic properties of the system over a wide range of temperature. 

Entropic sampling obtains when g(E) = D(E). This choice of g(E) renders P g {E) the 
same for all E, see Eq. (J7|). The system does a simple random walk on a one dimensional 
energy axis. All energy regions are visited with equal probability. The microstates on the 
paths that connect ordered and disordered phases in a first order phase transition would 
get equally sampled. A crucial issue that needs to be addressed pertains to the observation 
that we do not know D(E) a priori. We need a strategy to push g(E) closer and closer 
to D(E) iteratively. This we accomplish by employing Wang-Landau algorithm [ijj]. We 
divide the range of energy into a large number of equal width bins. We denote the discrete 
energy version of g(E) by the symbol {gi}- We start with gi = 1 V z; the subscript denotes 
energy bin index. We update {g,i} after every Monte Carlo step. Let us say the system 
visits a microstate in a Monte Carlo step and let the energy of the visited microstate fall 
in, say the m-th energy bin; then g m is updated to / x g m , where / is the Wang-Landau 
factor, see below. The updated {g^} becomes operative immediately for determining the 
acceptance/rejection criteria of the very next trial microstate. We set / = fo > 1 for the 
zeroth run. We generate a large number of microstates employing the dynamically evolving 
p. During the run we also build a histogram of energy of the miocrostates visited by the 
system. At the end of a run we check if the energy histogram has spanned the desired 
energy range and is reasonably flat. Flatter the histogram closer is {<?j} to {-Dj}, where 
{Di} is discrete energy representation of D(E). From one run to the next, the range of 
energy spanned by the system would increase. A run should be long enough to facilitate 
the system to span an energy range reasonably well and to render the histogram of energy, 
approximately flat. At the end of, say the v-th run, the Wang-Landau factor for the next 
run is set as f — /„ +1 = \fjl- After several runs, the Wang-Landau factor / would be very 
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close to unity; this implies that there would occur no significant change in {gi} during the 
run. For example with the square-root rule and /o = e, we have fio = exp(2~ 10 ) = 1.001. 
It is clear that / decreases monotonically with increase of the run index and asymptotically 
reaches unity. Wang and Landau have recommended the square-root rule; any other rule 
consistent with the above properties of monotonicity and asymptotic convergence to unity 
would do equally well. 

We take the output {gi}, which has converged reasonably to {A} (indicated by the 
uniformity of energy histogram) from the Wang-Landau Monte Carlo and carry out a single 
long entropic sampling run which generates microstates belonging to g ensemble. During 
this final production run we do not update {gi}- By un- weighting and re- weighting, see Eq. 
(IH]), of the microstates sampled from the ^-ensemble during the production run, we calculate 
the desired properties of the system as a function of (3. This is the strategy we shall follow 
in the Wang-Landau simulation of liquid crystal elastomer system, described in this paper. 

The usefulness of the Wang-Landau algorithm has been unambiguously demonstrated 
for system with discrete energy spectrum. However when we try to apply the technique to 
systems with continuous energy, there are serious difficulties that need careful considerations. 
In the present simulation we follow the modified Wang-Landau algorithm proposed in 13]. 
It essentially consists of treating an entire Wang-Landau run as a single iteration; {g{} 
obtained at the end of a Wang-Landau run is taken as an input for the next iteration. The 
convergence of {g{\ to {D{\ is monitored by the flatness of the histogram of energy measured 
in each iteration. Once we get a reasonably flat histogram we stop the iteration and start 
a production run with the final {g{} obtained in the last iteration. For more details of the 
simulation see (li 



IV. RESULTS AND DISCUSSIONS 



We consider an homogeneous lattice liquid crystal elastomer, with no external field and 
no stress: i.e. J it j = 1 V and \hi) = V i, a = and ji = 1. We consider a system with 
linear size L = 6. We take the initial microstate with all parallel to each other and 

A = Ao = 1. Let Co denote the initial microstate and Eq be its energy. We probe the system 
in an energy range —500 to 50 expressed in reduced unit. This energy range is divided into 
8000 equal width bins. Let v denote the energy bin index of the initial microstate. We set 
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gi = 1 + V i. In each Monte Carlo step a director is chosen at random and its orientation 
is changed to a trial orientation employing Barker's method It essentially consists 

of randomly selecting one of the three pre-fixed orthogonal axes and rotating the director 
about the chosen axis by a small angle A8 sampled randomly and uniformly between and 
0.02 radians. A trial strain parameter, X t is obtained from the current strain parameter 
Ao randomly following the prescription At = Ao + 0.01 x (£ — 0.5) where £ is a uniformly 
distributed independent random number between and 1. These two operations give us a 
trial microstate C t of the lattice elastomer. The acceptance of the trial microstate is based 
on entropic sampling described earlier. Once a microstate is selected we update g as per 
Wang-Landau algorithm, described earlier. We carry out in this fashion one Wang-Landau 



run; the g function at the end o 
run. The details are described in 



a Wang-Landau run is taken as an input for the next 



131 ] . Once we get an ' entropic function ' g that leads to 



approximately flat histogram of energy, we stop the iteration process and start a production 
run in which we obtain a large number microstates all belonging to g - ensemble. From the 
microstates sampled from the g ensemble we calculate all the desired macroscopic properties 
of the liquid crystal elastomer system, at various temperatures through un-weighting and 
re- weighting given by Eq. (Q. 

Figure (1) depicts orientational order S as a function of temperature T, for various values 
of coupling parameter 7 = 1, 0.8, 0.6, 0.4, and 0. For 7 = 1 we find that order parameter 
drops sharply when temperature increases. For smaller 7, the transition becomes less sharp 
and occurs at lower temperature. The loss of sharpness in the transition for small 7 is also 
due to the small system size considered in the study. Figure (2) depicts the strain parameter 
A as a function of temperature for various values of 7. At high temperature the system is 
isotropic and hence irrespective of the coupling strength, the strain is zero i. e A = 1 for all 
temperatures. When the temperature is lowered, say below T t , strain develops, see Fig. (2). 
The value of T t is higher for larger 7. For 7 = 1 the strain rises rather steeply with lowering 
of temperature. For lower values 7 the increase of strain with decrease of temperature for 
T < T t , is not large. For 7 = the orientational and elastic degrees of freedom are 
completely de-coupled and hence A = 1 for all temperatures, i.e. no strain develops even 
when temperature is lowered and nematic order sets in. These results are consistent with 
experimental observations. 

Figure (3) depicts specific heat, calculated from the fluctuations of energy, as a function 
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of temperature for various values of 7. For 7 = 1, we find that the specific heat shows a 
sharp maximum at the transition temperature. For lower values of 7 the transition occurs 
at lower temperatures. This is in agreement with the mean-field arguments and Uchida [20] ; 
also the peak is broader at lower 7. 

Figure (4) depicts the transition temperature as a function of 7. The transition tem- 
perature is taken as the value of T at which the specific heat shows a maximum. As 7 
increases the transition temperature increases slowly initially; when 7 increases beyond 0.6 
the transition temperature increases rather steeply. 

As seen from Figs. (1) and (2), both S and A increase with decrease of temperature. To 
see the nature of their correlation we have plotted in Fig. (5), A versus S for various values 
of 7. For 7 = the strain is zero (A = 1) and is independent of S as expected since in this 
case the orientational and elastic degrees of freedom are uncorrelated. For small values of 7 
the strain parameter A scales linearly with S over the full range of temperature. For 7 = 1, 
the scaling is linear for S < 0.6 and A < 1.25 which correspond to t emp erature greater than 
about 0.8. This is consistent with the results of earlier simulation 11[ which showed linear 
scaling between A and S. However, we find that for lower temperatures the scaling of A 
with S is nonlinear. The strain increases steeply to large values as the orientational order 
increases and attains its maximum value of unity. 

One of the advantages of multicanonical Monte Carlo techniques is that we can calculate 
density of states up to a normalization. The array {gi} would give an approximate estimate of 
density of states if the corresponding histogram of energy of sampled microstates is relatively 
flat. Fig. (6) depicts g function and the corresponding energy histogram. The energy 
histogram is relatively flat. From the g function we can calculate the microcanonical entropy 
and free energy up to an additive constant. Note that entropy or free energy calculations 
are very difficult if not impossible, from conventional Metropolis Monte Carlo simulations. 
We show in Fig. (7) variation of (relative) free energy with energy for 7 = 1 at three 
temperatures bracketing the transition point. We see clearly that the transition is first 
order and strong. Figure (8) depicts free energy versus E for 7 = 0.8. The transition is still 
first order but is relatively weak. For smaller values of 7 the transition weakens further. For 
7 less than 0.4 or so, the free energy barrier is not discernible. 
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V. CONCLUSIONS 



We have reported in this paper results on phase transition in liquid crystal elastomers 
employing multicanonical Monte Carlo simulations. We have considered the nature of 
transition for various strengths of coupling of the local nematic director to a global strain. 
We find that at maximal coupling, the transition is strongly first order. As the coupling 
becomes weaker, the transition becomes weakly first order. We find that the scaling of 
nematic order with strain is non-linear especially when the local nematic directors are 
strongly coupled to global strain and at low temperatures. 
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FIG. 1: Orientational ordre parameter S versus temperature for various values of the the coupling 
parameter 7. 
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FIG. 2: Strain parameter A versus temperature for various values of the the coupling parameter 7. 
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FIG. 3: Specific heat Cy versus temperature for various values of the coupling parameter 7; Cy 
has been obtained from the energy fluctuations. 
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FIG. 4: Transition temperature versus the coupling parameter 7; the temperature at which the 
specific heat exhibits peak is taken as the transition point. 
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FIG. 6: (Top) The g function that approximates the density of states D; logarithm of g gives the 
microcanonical entropy; what is plotted is logarithm of entropy as a function of ewnergy. (Bottom) 
The histogram of energy of mictrostates sampled during the production run employing entroping 
sampling with acceptance determined by the density of states depicted in the top. 
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FIG. 7: Free energy versus energy for temperatures above, below and at the transition point for 
the system with 7 = 1.0. It is clear that the transition is first order. 
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